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ABSTRACT 

We present results of a Bayesian analysis of radial velocity (RV) data for the star HIP 5158, 
confirming the presence of two companions and also constraining their orbital parameters. 
Assuming Keplerian orbits, the two-companion model is found to be e^^ times more probable 
than the one-planet model, although the orbital parameters of the second companion are only 
weakly constrained. The derived orbital periods are 345.6 ± 2.0 d and 9017.8 ± 3180.7 d 
respectively, and the corresponding eccentricities are 0.54 ± 0.04 and 0.14 ± 0.10. The limits 
on planetary mass (msini) and semimajor axis are (1.44 ± 0.14 A/j, 0.89 ± 0.01 AU) and 
(15.04 ± 10.55 Afj, 7.70 ± 1.88 AU) respectively. Owing to the large uncertainty on the mass 
of the second companion, we are unable to determine whether it is a planet or a brown dwarf. 
The remaining 'noise' (stellar jitter) unaccounted for by the model is 2.28 ± 0.31 m/s. We 
also analysed a three-companion model, but found it to be times less probable than the 
two-companion model. 
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1 INTRODUCTION 

Extrasolar planetary research has made great advances in the last 
decade as a result of the data gathered by several ground and space 
based telescopes and so far more than 500 extrasolar planets have 
been discovered. More and more planets with large orbital peri- 
ods and small velocity amplitudes are now being detected due to 
remarkable improvements in the accuracy of RV measurements. 
With the flood of new data, more powerful statistical techniques 
are being developed and applied to extract as much information 
as possible. Traditionally, the orbital parameters of the planets and 
their uncertainties have been obtained by a two stage process. First 
the period of the planets is determined by searching f or periodicity 
in the RV dat a using the Lomb-Scargle periodogram dLombll 19761 : 
IScarglelll982h . Other orbital parameters are then determined using 
minimisation algorithms, with the orbital period of the planets fixed 
to the values detennined by Lomb-Scargle periodogram. 

Bayesian methods have several advantages over traditional 
methods, for example when the data do not cover a complete or- 
bital phase of the planet. Bayesian inference also provides a rig- 
orous way of performing model selection which is required to de- 
cide the number of planets favoured by the data. The main prob- 
lem in applying such Bayesian model selection techniques is the 
computational cost involved in calculating the Bayesian evidence. 
Nonetheless, Bayesian model selection has the potential to improve 
the interpretation of existing observational data and possibly detect 
yet undiscovered planets. Recent advances in Marko-Chain Monte 
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Carlo (MCMC) techniques (see e.g. lMackavl[2003h have made it 
possible for Bayesian techniques to b e applied to extrasolar plane- 
tary searches (se e e.g. 'Gregorv"2005^ 'Ford"2005l iFord & Gregorvl 
l2007l : lBalan&L ahav 2009). Feroz et al. (2010) recently presented 
a new Bayesian method for determining the number of extrasolar 
planets, as well as for inferring their orbital parameters, without 
having to calculate directly the Bayesian evidence for models con- 
taining a large number of planets, although this method is not re- 
quired for the analysis of the simple system considered here. 

HIP 5158 is a K-type main sequence star at a distance o f 
45 ± 4 pc with mass 0.780 ± 0.021 Mq fco Curto et Zll2010h . 
In this paper, we present a Bayesian analysis of the existing, high- 
precis ion RV measureme nts of HI P 5158 from the HAR PS instru- 
ment jMavor et all2003h . given in lLo Curto et all (|2010'), who de- 
tected a planet with period 345.72 ± 5.37 d. They did not find any 
correlation of the RV with the bispector span of the cross correla- 
tion function for HIP 5158. They also did not find any periodicity 
in the stellar activity indicator log -RjjK' making it unlikely for any 
long period RV variability to be caused by stellar activity. More- 
over, they noticed an additional quadratic drift in the RV data which 
they inferred to indicate the presence of another body orbiting the 
star, but they were not able to constrain its orbital parameters. 

The outline of this paper is as follows. We give a brief intro- 
duction to Bayesian inference in Sec.|2]and describe our method for 
calculating the number of planets favoured by the data. Our model 
for calculating RV is described in Sec.|3] In Sec.|4]we describe our 
Bayesian analysis methodology. We apply our method to RV data 
sets on HIP 5158 in Sec.[5]and present our conclusions in Sec.|6l 
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2 BAYESIAN INFERENCE 

Our planet finding methodology is built upon the principles of 
Bayesian inference, and so we begin by giving a brief sum mary 
of this framework. We refer the reader to iFeroz et all fcoiol) for a 
more thorough discussion on Bayesian object detection methods. 

Bayesian inference methods provide a consistent approach to 
the estimation of a set of parameters in a model (or hypothesis) 
H for the data D. Bayes' theorem states that 



Pr(0|D,ii') = 



Pr(D| &,H)Pr{@\H) 



(1) 



Pr(D|/f) 

where Pr(0jD,_ff) = -P(0) is the posterior probability distri- 
bution of the parameters, Pr(Dj©, H) = £{&) is the likelihood, 
Pi{@\H) = 7r(0) is the prior, and Pr(D|//) = 2 is the Bayesian 
evidence. 

In parameter estimation, the normalising evidence factor can 
be ignored, since it is independent of the parameters 0, and in- 
ferences are obtained by taking samples from the (unnormalised) 
posterior using standard MCMC methods, where at equilibrium 
the chain contains a set of samples from the parameter space dis- 
tributed according to the posterior This posterior constitutes the 
complete Bayesian inference of the parameter values, and can be 
marginalised over each parameter to obtain individual parameter 
constraints. 

In contrast to parameter estimation problems, for model se- 
lection the evidence needs to be evaluated and is simply the factor 
required to normalize the posterior over 0: 



(2) 



where D is the dimensionality of the parameter space. As the av- 
erage of the likelihood over the prior, the evidence is larger for 
a model if more of its parameter space is likely and smaller for a 
model with large areas in its parameter space having low likelihood 
values, even if the likelihood function is very highly peaked. Thus, 
the evidence automatically implements Occam's razor; a simpler 
theory with compact parameter space will have a larger evidence 
than a more complicated one, unless the latter is significantly bet- 
ter at explaining the data. The question of model selection between 
two models Hq and Hi can then be decided by comparing their 
respective posterior probabilities given the observed data set D, as 
follows 



R 



Pr(J/ijD) _ Pr(Dj//i)Pr(iyi) _ Zi Pr{Hi) 



Pr{Ho\T)) Pr(Dj//o)Pr(i?o) ZoPr(Ho) 



(3) 



where Pr{Hi)/ Pt{Ho) is the a priori probability ratio for the two 
models, which can often be set to unity but occasionally requires 
further consideration. The natural logarithm of the ratio of poste- 
rior model probabilities (sometimes termed the posterior odds ra- 
tio) provides a useful guide to what constitutes a significant differ- 
ence between two models: 



Alnii = In 



Pr(/fi|D) 



Pr(//o|D) 



= In 



Zi Pr(Hi) 



Zo Pr{Ho) 



(4) 



The evaluation of the number of objects favoured by the data 
is a model selection problem. By considering a series of models 
/fjVobj ' ^^ch with a fixed number of objects, i.e. with A^obj ~ 
0, 1, 2, . . ., one can infer A'^obs by identifying the model with the 
largest marginal posterior probability Pr(//jv^jj. |D). The proba- 
bility associated with A'^nbj = is often called the 'null evidence' 
and provides a baseline for comparison of different models. Indeed, 
this approach has been adopted previously in exoplanet studies (see 



e.g. lGregorv & Fischeij fcoiol) ). albeit using only lower-bound es- 
timates of the Bayesian evidence for each model. 

Evaluation of the multidimensional integral in Eq.|2]is a chal- 
lenging numerical task. Standard techniques like thermodynamic 
integration are extremely computationally expensive which makes 
evidence evaluation at least an order of magnitude more costly than 
pa rameter estimat ion. The nested sampling approach, introduced 
bv lSkillind l2004h. is a Monte Carlo method targeted at the effi- 
cient calculation of the e vidence, but also produ ces posterior in- 
ferences as a by-product. iFeroz & Hob son (2008) and iFeroz et al.l 
( 2009b) built on this nested sampling framework and have recently 
introduced the MultiNest algorithm which is very efficient in 
sampling from posteriors that may contain multiple modes and/or 
large (curving) degeneracies and also calculates the evidence. This 
technique has greatly reduces the computational cost of Bayesian 
parameter estimation and model selection and has already been ap- 
plied to several model se lections problem in astrophysics (see e.g. 
iFeroz et al .l2008ll2009cllah . We employ this technique in this paper. 



3 MODELLING RADIAL VELOCITIES 

Observing planets at interstellar distances directly is extremely dif- 
ficult, since the planets only reflect the light incident on them from 
their host star and are consequently many times fainter. Nonethe- 
less, the gravitational force between the planets and their host star 
results in the planets and star revolving around their common cen- 
tre of mass. This produces doppler shifts in the spectrum of the 
host star according to its RV, the velocity along the line-of-sight to 
the observer. Several such measurements, usually over an extended 
period of time, can then be used to dete ct extrasolar planets. 

Following the formalism given in iBalan & Lahavl ( l2009h . for 
TVp planets and ignoring the planet-planet interactions, the RV at 
an instant ti observed at jth observatory can be calculated as: 

iVp 

■"{tij) = Vj-'^Kp [sin(/i,p + zup) + Cp sin(ti7p)] , (5) 
p=i 



where 

Kp 



systematic velocity with reference to jth observatory, 
velocity semi-amplitude of the pth planet, 



ujp = longitude of periastron of the pth planet, 

/i,p — true anomaly of the pth planet, 

Bp = orbital eccentricity of the pth planet, 

start of data taking, at which periastron occurred. 

Note that /i,p is itself a function of gp, the orbital period Pp of the 
pth planet, and the fraction Xp of an orbit of the pth planet, prior to 
the start of data taking, at which periastron occurred. While there 
is a unique mean line-of-sight velocity of the center of motion, it is 
important to have a different velocity reference Vj for each obser- 
vatory/spectrograph pair, since the velocities are measured differ- 
entially relative to a reference frame specific to each observatory. 

We also model the intrinsic stellar variability s ('jitter'), as a 
source of uncorrelated Gaussian noise in addition to the measure- 
ment uncertainties. Therefore for each planet we have five free pa- 
rameters: K, tjj, e, P and x- In addition to these parameters there 
are two nuisance parameters V and s, common to all the planets. 

These orbital parameters can then be used along with the stel- 
lar mass rris to calculate the length a of the semi-major axis of the 
planet's orbit around the centre of mass and the planetary mass m 
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Pai'ameter 


Prior 


Matliematical Form 


Lower Bound 


Upper Bound 


P (days) 


Jeffreys 


1 

Pln(Pmax/P„i„) 


0.2 


365, 000 


K{m/s) 


Mod. Jeffreys 





-ffmax(Pmi„/Pi)'''3(l/yi - e?) 


ln(l+ (K,„ax /if o) (-PminZ-Pi ) ^ (1/ \/l - ) ) 


V (m/s) 


Uniform 


1 

Vmin-Vmax 




-^max 


e 


Uniform 


1 





1 


zu (rad) 


Uniform 


1 

27r 





27r 


X 


Uniform 


1 





1 


s (m/s) 


Mod. Jeffreys 


(s + so)"^ 
ln(l + s„,ax/so) 





^max 



Table 1. Prior probability distributions. 



as follows: 



Os sm I 



KPVT 



2-K 



m sm t ~ 5 , 

(27rG)3 

msassini 

a « (8) 

m sm t 

where as is the semi-major axis of the stellar orbit about the centre- 
of-mass and i is the angle between the direction normal to the 
planet's orbital plane and the observer's line of sight. Since i can- 
not be measured with RV data, only a lower bound on the planetary 
mass m can be estimated. 



4 BAYESIAN ANALYSIS OF RADIAL VELOCITY 
MEASUREMENTS 

There are several RV search programmes looking for extrasolar 
planets. The RV measurements consist of the time ti of the ith 
observation, the measured RV Vi relative to a reference frame and 
the corresponding measurement uncertainty cti . These RV measure- 
ments can be analysed using Bayes' theorem given in Eq.[T]to ob- 
tain the posterior probability distributions of the model parameters 
discussed in the previous section. We now describe the form of the 
likelihood and prior probability distributions. 

4.1 Likelihood function 

As discussed in iGregorvl JioOTh , the eiTors on RV measurements 
can be treated as Gaussian and therefore the likelihood function 
can be written as: 

{v{e;ti)~viy 





/Vp 


In 2 


s (m/s) 


(6) 


1 


30.05 ±0.14 


10.31 ± 1.09 




2 


78.19 ±0.15 


2.28 ± 0.31 


(7) 


3 


70.68 ±0.16 


2.28 ± 0.28 



£(9) = n 



: exp 



(9) 



where Vi and ai are the i*" ^ RV measurement and its corresponding 
uncertainty respectively, v{0;ti) is the predicted RV for the set of 
parameters Q, and s is intrinsic stellar variability. A large value 
of s can also indicate the presence of additional planets, e.g. if a 
two-planet system is analysed with a single-planet model then the 
velocity variations introduced by the second planet would act like 
an additional noise term and therefore contribute to s. 



4.2 Ciioice of priors 

For parameter estimation, priors become largely irrelevant once the 
data are sufficiently constraining, but for model selection the prior 



Table 2. The evidence and jitter values for the system HIP 5158. The null 
evidence (—255.3) has been subtracted from each In 2 value. 



dependence always remains. Therefore, it is important that priors 
are selected base d on physical con siderations. We follow the choice 
of priors given in lOregorvl ( l2007h . as shown in Table [T] 
The modified Jeffreys prior. 



Pr{e\H) = 



(10) 



behaves like a uniform prior for 9 <C 9o and like a Jeffreys prior 
(uniform in log) for 2> ^o- We set Ko = so = 1 ni/s and 
Jfmax ~ 2129 m/s, which corresponds to a maximum planet-star 
mass ratio of 0.01. 



5 RESULTS 

As discussed in Sec.|2l the number of companions orbiting a star 
can be determined by analysing the RV data, starting with A^p — 
and increasing it until the Bayesian evidence for A'p = n com- 
panions starts to drop off. The resulting evidence and jitter values 
for HIP 5158 RV data, after subtracting a mean RV of 14.9105 
km/s from it, are presented in Table |2] We can clearly see A'^p — 2 
is the strongly favoured model. Adopting the 2-companion model, 
the estimated parameter values are listed in Table |3] while the 1-D 
marginalised posterior probability distributions are shown in Fig.[T] 
As mentioned in Sec.[T] stellar activity is unlikely to be the source 
of any long period RV variability of HIP 5158. Moreover, the in- 
ferred period of HIP 5158 c is more than two orders of magnitude 
hi gher than the rotation period of the star (Prot(log -Rhk) = 42.3 
d. Ilo Curto et al.ll2OI0l) . We therefore conclude that it is unlikely 
for the second signal to be generated by the photospheric activity 
of HIP 5158. 

The mean RV curve for 1 and 2-companion models are over- 
laid on the RV measurements in Fig.|2] The RV residuals after fit- 
ting for 1 and 2 companions along with nuisance parameters (con- 
stant systematic velocity and stellar jitter) are shown in Fig.|3] and 
are obtained by subtracting the estimated mean RV from the mea- 
sured RV and adding the estimated stellar jitter in quadrature to the 
quoted RV uncertainties. The 1 -companion fit is clearly very poor 
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Figure 1. 1-D marginalised posterior probability distributions for the parameters of the two bodies found orbiting HIP 5158. 



due to the presence of quadratic drift in the RV data, whereas the 
residuals from the 2-companion fit are consistent with a noise-only 
model. Fig. |4] shows histograms of normalised data residuals, de- 
fined as follows: 



V, — V\ 

where iii, ai and vi are the measured RV, measured RV uncer- 
tainty and the estimated mean RV respectively, all at time ti and 
s is the estimated mean stellar jitter. The histogram for the nor- 
malised residuals in the 1 -companion fit (left panel) has a similar 
width to the corresponding histogram for the 2-companion fit (right 
panel), but it should be noted that s appearing in the normalisation 
of Eq.[TT]is much smaller in the latter case (see Table |2) 

Comparin g our parameter values with those given in 
ILo Curto et all feoiot) . we see that our parameter estimates for 
planet HIP 5158 b are in very good agreement. The large uncer- 
tainties on the parameter estimates of HIP 5158 c are because not 
enough phase has been covered by the RV measurements. Never- 
theless, we can still determine that the period of the second com- 
panion is greater than 10 years with the 95% confidence interval 
being (3567, 18194) days. The corresponding 95% confidence in- 
tervals for minimum mass and semi-major axis of the second com- 
panion are (2.3,41.5) A/j and (4.2,12.5) AU respectively and 
therefore we can not determine whether it is a planet or a brown 
dwarf. 



Parameter 


HIP 5158 b 


HIP 5158 c 


P (days) 


345.63 ± 1.99 


9017.76 ± 3180.74 




(341.82,349.53) 


(3567.13, 18193.94) 




(344.74) 


(11135.63) 


A" (m/s) 


59.05 ± 7.73 


170.54 ± 110.17 




(49.89, 77.47) 


(36.59,416.85) 




(56.65) 


(218.49) 


e 


0.54 ±0.04 


0.14 ±0.10 




(0.47, 0.63) 


(0.01,0.40) 




(0.54) 


(0.00) 


■cu (rad) 


1.23 ± 0.07 


2.48 ± 1.32 




(1.06, 1.34) 


(0.33,5.90) 




(1.24) 


(2.23) 


X 


0.02 ± 0.00 


0.57 ±0.23 




(0.02,0.02) 


(0.12,0.95) 




(0.02) 


(0.59) 


msini (Mj) 


1.44 ±0.14 


15.04 ±10.55 




(1.26, 1.77) 


(2.34,41.49) 




(1.34) 


(19.53) 


a(AU) 


0.89 ± 0.01 


7.70 ± 1.88 




(0.87,0.91) 


(4.22, 12.50) 




(0.87) 


(8.83) 



Table 3. Estimated parameter values for the two bodies found orbiting HIP 
5158. The estimated values in the first row are quoted as /i ± cr where 
fi and a are the posterior mean and standard deviation respectively. The 
interval in the second row is the 95% confidence interval and the number in 
parenthesis in the third row is the maximum-likelihood parameter value. 
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Figure 2. RV measurements, with Icr eiTorbars, and the mean fitted RV 
curve with 1 (red solid line) and 2-companions (green dashed line) for HIP 
5158. 
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Figure 3. RV residuals, with Icr errorbars that include the estimated stellar 
jitter, after fitting for one (red sohd bars) and two companions (blue dotted 
bars). 



6 CONCLUSIONS 

We have presented a Bayesian analysis of HIP 5158 RV data. Us- 
ing Bayesian model selection, we have found a very strong signal 
for the presence of two companions orbiting HIP 5158. Our esti- 
mated orbital parameters for the pla net HIP 5158 b are in excellent 
agreement with the values given in ILo Curto et all ( I2OI0I) . We de- 
termined the orbital period of HIP 5158 c to be at least 10 years but 
the presence of large uncertainties in its orbital parameters do not 
allow us to determine whether it is a planet or a brown dwarf. 
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Figure 4. Histograms of normalised RV residuals, as defined in Eg. II II 
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